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Abstract 

This paper presents a survey of a set of simulations of accretion disks orbiting a 
rapidly rotating Kerr black hole and embedded in a large-scale initial magnetic field. 
Each simulation uses a common state for the initial torus, including an MRI seed-field 
consisting of poloidal loops along isodensity contours, and differ only in the strength 
of the initial large-scale magnetic field in relation to these poloidal loops. Simulations 
with a weak large-scale initial field differ little from simulations that evolve from only 
the poloidal field. Simulations where the large-scale initial field distorts or completely 
overwhelms the poloidal loops show more extensive regions of turbulence, due to the 
action of the MRI on the large-scale field. However, the overall structure of the late- 
time state of the simulations remains qualitatively unchanged, that is each simulation 
sees the emergence of a turbulent accretion disk, axial jets, a funnel-wall outflow and 
an extensive corona. 



1 Introduction 



The numerical study of black hole accretion has been an area of ongoing research, 
beginning with the pioneering studies of Wilson (1972). Simulations of accretion driven 
by the Magneto-Rotational Instability (MRI; Balbus & Hawley, 1998) also have a long 
history, with various approximations used to mimic the properties of the central mass 
(e.g., the pseudo-Newtonian simulations of Hawley & Krolik, 2001). More recently, 
codes such as the General Relativistic MHD (GRMHD) code of De Villiers & Hawley 
(2003; DH03a) have been used to study accretion disks orbiting Kerr black holes and 
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the jets that arise self-consistently from the accretion flow (De Villiers, Hawley & 
Krolik, 2003; DHK). The key result from these simulations is that GRMHD is required 
to correctly capture the effect of black hole spin on the energy of the jets (De Villiers 
et al., 2005; D05a), which are powered by an MHD interaction between the accretion 
flow and the spinning black hole. 

Another line of inquiry deals with simulations of radial accretion onto black holes. 
In such simulations, for instance those described by Komissarov (2005), an initial state 
consisting of a tenuous plasma embedded in a large-scale magnetic field is evolved. 
Such simulations feature many interesting aspects, of immediate interest are the wind- 
ing of the magnetic field by black hole rotation, and the generation of broad, energetic 
outflows. Such simulations seem to encounter a generic problem: the outflow which de- 
velops tends to be uncollimated and severely limits further accretion since any infalling 
gas at large radii is soon blown out by the outflow produced near the black hole. 

Accretion disk simulations do not suffer from this drawback: the initial torus in the 
equatorial plane serves as a long-lived mass reservoir that can sustain axial outflows for 
very long periods (D05a). However, most of the accretion disk simulations carried out 
with the GRMHD code have not featured a large-scale initial magnetic field, and instead 
trigger the MRI in an initial torus by using a weak initial magnetic field confined to 
the region immediately surrounding the pressure maximum of the torus. Since ambient 
fields are likely to exist in some form in an astrophysical context, it is not simply of 
academic interest to consider how the evolution of an accretion disk would be affected 
if the initial torus were immersed in a large-scale weak magnetic field. In the context of 
GRMHD simulations, this is readily achieved using analytic solutions for a large-scale 
field due to Wald (1974; W74). 

Attempts at merging the two approaches with a general relativistic code have been 
reported. Nishikawa et al. (2005) evolve a thin accretion disk in the Schwarzschild 
metric in three spatial dimensions. However, those simulations use an extremely strong 
initial magnetic field (ft of order unity in the initial disk), which may pose problems 
for the growth of the MRI (Hawley & Balbus, 1992; HB92), and are of very short 
duration, on the order of a few percent of the typical duration reported here, again 
precluding proper development of the MRI. De Villiers et al. (2005; D05b) reported 
three axisymmetric simulations and one 3D simulation of an accretion disk orbiting 
a Kerr black hole evolved from an initial torus immersed in a weak large-scale field; 
the i? v f simulations, though applied to the collapsar scenario, are analogous to the 
weak-field simulation (the W simulation; see below) discussed here. 

This paper presents a family of axisymmetric (2.5D) simulations which evolve an 
initial torus immersed in a Wald magnetic field and orbiting a Kerr black hole. The 
simulations differ only in the relative strength of the large-scale field that is superim- 
posed on the traditional MRI seed field, that of poloidal loops embedded in the initial 
torus. Technical details are left to the appendices. 
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2 Overview of the Simulations 



To simulate accretion disks, begin with an initial equilibrium torus in the space-time 
of a rotating black hole. This equilibrium torus is a differentially rotating system, with 
a near-Keplerian distribution of angular momentum. In order for this torus to evolve 
into an accretion disk, it needs to be perturbed from its equilibrium state by intro- 
ducing small density fluctuations in the torus, and, most importantly, by supplying 
a weak magnetic field that constitutes the seed field which will trigger the Balbus- 
Hawley magneto-rotational instability. The MRI amplifies the seed magnetic field by 
shearing. The MRI is characterized by a linear growth phase of the low-order modes 
of the magnetic field, a process which take place on the time scale of the orbital pe- 
riod of the fluid. This is followed by an exponential growth phase due to non-linear 
interaction between the growing modes. This non- linear interaction gives rise to the 
turbulent transport of angular momentum, transforming the initial torus into an accre- 
tion disk. In axisymmetry, the process of turbulent transport of angular momentum is 
not self-sustaining, due to Cowling's anti-dynamo theorem; axisymmetric simulations 
have shown that turbulence dies down, typically within a few orbits of onset. 

In earlier 3D simulations (e.g. DHK), the MRI seed field was confined to the region 
near the pressure maximum of the initial torus. Two configurations have been used: 
poloidal loops laid down along isodensity contours and, less frequently, large-scale 
toroidal loops. In such simulations, the MRI evolves on the local time-scale dictated 
by the orbital period at the initial pressure maximum, since this is the early locus of 
the MRI. Therefore, the orbital period at the initial pressure maximum constitutes 
a natural measure of simulation time for these simulations. In the 2.5D simulations 
discussed here, the MRI seed field consists of poloidal loops on which are superimposed 
a pervasive large-scale field whose overall orientation is parallel to the spin-axis of 
the black hole. Since this large-scale field threads the entire initial torus, the orbital 
period at the pressure maximum is not used as a time measure; comparisons between 
simulations will be done using the common relativistic time unit of black hole mass. 

Four simulations are discussed in this paper. The simulations differ only in the 
strength and orientation (parallel or anti-parallel to the spin axis) of the large-scale 
initial magnetic field. Each simulation features a torus with structural parameters 
listed in Table ^ The initial tori are immersed in a tenuous medium which consists of 
a cold, radially infalling dust. The initial pressure maximum of these tori is at a radius 
or r/M = 16.1 in the equatorial plane, which is relatively close to the black hole; the 
main reason for this choice is that it allows rapid evolution of the accretion disk. The 
initial tori are perturbed by introducing random density fluctuations (1% amplitude) 
as well as a weak magnetic field in the form of poloidal loops with an average (5 (the 
ratio of gas to magnetic pressure) of 100. The initial magnetic field is then completed 
by superimposing a large-scale magnetic field on the poloidal loops; the strength of 
the large-scale field is set by a scaling constant, Bq in Wald's solution, and differs for 
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each model. The strength of the large-scale field is measured in relation to the initial 
poloidal loops: 

• Weak field: where the large-scale field is much weaker in the torus than the 
poloidal loops and barely perturbs the structure of the loops (this simulation is 
denoted W); 

• Moderate field: where the large-scale field is comparable in magnitude through 
the torus to the poloidal loops; two cases are considered, where the field is oriented 
in the "up" direction (M u ), and the "down" direction (Md), and depending on the 
orientation, this field cancels the "vertical" portion of the poloidal loops either 
inward or outward of the pressure maximum; 

• Strong field: where the large-scale field dominates the poloidal loops, effectively 
washing out these loops in the torus (denoted S). The term "strong field" is a 
relative one; the strength of this field was adjusted to just wash out the initial 
poloidal loops, and the overall strength of the field, as measured by the plasma 
/^-parameter, (the ratio of gas to magnetic pressure) remains weak. 

Each of these initial states is shown in Figure Q where gas density is shown along with 
arrows showing the orientation of the magnetic field. The effect of the large-scale field 
on the poloidal loops can be readily seen in these simple diagrams. In more quantitative 
terms, typical values of (5 in the initial torus lie in the range of 100 to 1000, for all 
models: the initial field in the region where the MRI will operate is weak. The values 
of j3 in the region filled with a tenuous radial dust range from 10 -3 to 10~ 5 , with the 
lower values found in the S model. 

Each simulation is run for a duration of at least 2000 M. This duration is sufficient 
to pass through the initial transient phase where initial torus becomes accretion disk 
(during the first 600 M or so of simulation time), and resolve several complete orbits of 
the accretion disk, at which time the turbulent transport of angular momentum begins 
to die down (in axisymmetry) . Complete dumps of the code variables are saved every 
2 M increment in simulation time. Shell-integrated diagnostic values are computed and 
saved every 1 M increment in simulation time. 
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Table 1: Torus Simulation Parameters. 



Black hole: 






spin 


a/M 


0.90 


horizon 


r h /M 


1.44 


marginally stable orbit 


r ms /M 


2.32 


Torus: (See Note a) 






inner edge 


r in /M 


9.5 


pressure maximum 




16.1 


outer edge 


Tont/M 


36.0 


maximum height 


H/M 


9.0 


orbital period (at rp max ) 


T orh /M 


422 


T™\ J T~> 1 "I / r . 717-|7\ 

Dust Background: (See Note b) 






density contrast 


Pdust/ Pdisk 


10~ 6 


energy contrast 


Cdust/edisk 


io- n 


Grid: (See Note c) 






Zones 


N r x N e 


192 x 192 


Inner boundary 


r mi jM 


1.45 


Outer boundary 


'"max / M- 


120 


Polar axis offset 


A9/n 


lO" 3 


External Field: (See Appendix^ 






"Zero" field (Z) 


Bo 





Weak field (W) 


B 


+io- 8 


Moderate field (M u ) 


Bo 


+io- 2 


Moderate field (M d ) 


Bo 


-lo- 2 


Strong field (S) 


Bo 


+10+ 2 



Note a: Equations given in Appendix 151 K — 0.01, q = 1.68, and T = 4/3. Torus seeded with poloidal 
magnetic field loops with (3 = (-P g as)/(-fmag) ~ 100. Torus also given small (1%) density perturbation. 
Note b: Equations in Hawlcy, Smarr & Wilson (1984); scaled to grid-averaged density and energy contrast. 
Note c: Radial grid: cosh-scaling, outflow boundaries; #-grid: exponential scaling, reflecting boundaries. 
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Figure 1: Initial state for the four simulations. The initial torus is plotted using the logarithm of 
density (grey scale). Arrows show the local direction of the initial magnetic field, which is purely 
poloidal. The initial states are, clockwise from bottow left: W, M u , M^, S. 
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3 Summary of Numerical Evolution 



This section outlines the role of the large-scale initial magnetic field in the simulations. 
To better appreciate the interaction between the large-scale field and the initial poloidal 
loops, a brief digression is in order to recapitulate the simulations where each type of 
field is considered separately. 

In simulations where the MRI is triggered solely by initial poloidal loops, the trans- 
formation of initial torus into accretion disk proceeds along the lines discussed exten- 
sively by DHK: the initial rapid growth of the toroidal magnetic field by shear during 
the first orbit of the initial torus; the saturation of the poloidal field MRI, which 
progresses outward through the torus; and the evolution by MHD turbulence of the 
accretion disk. The late-time structure can be divided into five regions: 

• The main body of the disk, a turbulent, massive wedge. Gas pressure dom- 
inates here and the magnetic and velocity fields are highly tangled. The outer 
part of the disk moves radially outward with time as it gains angular momentum. 

• The coronal envelope, a region of low density surrounding the disk, where gas 
and magnetic pressure are comparable. 

• The inner torus and plunging region. The inner torus, lying just outside 
the marginally stable orbit, temporarily accumulates gas accreting from the main 
disk. The plunging region, where accreting gas spirals into the black hole, lies 
inward of the inner torus. 

• The funnel- wall jet, an outflow along the centrifugal barrier which originates 
near the inner torus. The density here is small compared to the disk, but one 
to two orders of magnitude greater than in the axial funnel. The intensity of 
the jet is spin-dependent: the funnel-wall jet is strongest when black hole spin is 
near-extremal. The jet is also fed by entrainment from the adjacent corona. 

• The axial funnel, a magnetically-dominated region with tenuous, hot, relativistic 
outflow. The magnetic field is predominantly radial. 

Perhaps the most important point in the present context is that the magnetic field 
in these simulations, though initially confined to the torus, eventually populates the 
entire simulation volume (Hirose et al, 2004); the MRI-driven accretion flow produces 
its own large-scale field. DHK and companion papers documented full 3D simulations; 
simulations carried out in axisymmetry tend to produce a more violent transient phase 
and the turbulent disk exhibits the channel solution of HB92. However, in both 2.5D 
and 3D simulations, the above facts remain essentially unchanged with one important 
difference: in axisymmetry, the MRI is not self-sustaining and the rate of accretion 
dies down within about ten orbital periods (at the initial pressure maximum). 
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In simulations where a large-scale magnetic field interacts with a radially infalling 
dust, the process which drives the observed unbound outflows is the advection and 
subsequent winding up of the initial magnetic field through frame dragging. The time- 
scale for this process is very short compared to the accretion disk simulations: for a 
Kerr black hole with a spin of a/M = 0.9, as is used in the present simulations, the 
process of turning the initial infall into a broad, unbound outflow takes place within 
approximately 100 M of simulation time. Within this short period, the density of the 
gas in the vicinity of the black hole drops precipitously. The simulation results become 
unreliable when the density and energy approach the numerical floor of the code, at 
approximately 300 to 400 M in simulation time for the GRMHD code (De Villiers, 
unpublished simulations). As mentioned earlier, the lack of a stable mass reservoir 
near the black hole dooms such simulations to a short life. 

3.1 The Weak-field Case 

The evolution of the W simulation is, in all respects, very similar to a benchmark 
simulation with no initial large-scale field (referred to from hereon as the Z, or zero 
large-scale field simulation) . Figure [21 presents the late-time density distribution as 
well as the orientation of the poloidal component of the magnetic field for the W and 
Z simulations. The overall size of the main disk, the shape of the corona, the presence 
of an inner torus just outside the marginally stable orbit, and the poloidal magnetic 
field structure are virtually identical. As will be seen later, other diagnostics are also 
in close agreement for these two simulations. 




Figure 2: Late-time structure (t/M = 1600) of Weak-field, W, simulation (left panel) compared 
to the "zero" field, Z, simulation (right panel). Gas density is plotted on a common logarithmic 
scale. The arrows show the local direction of the poloidal component of the magnetic field. Each 
panel spans a physical extent of 50 M; the black hole is centered on the left edge of each panel. 
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3.2 The Moderate-field Cases 



Differences with the case with no initial large-scale field begin to emerge in the M u 
and Md simulations. As mentioned above, the time-scale for the winding up of the 
large-scale field is very short; in simulations with only this field and radial dust, an 
unbound outflow develops promptly near the black hole and impedes further accretion. 
In the moderate- and strong-field simulations presented here, this behaviour persists, 
but coexists with the early growth phase of the MRI in the torus, which extrudes a thin, 
dense current sheet along the equatorial plane. The clash between these two events is 
captured in Figure EJ which shows the "shoving match" between the nascent accretion 
flow and the transient outflow from the large-scale field. A sharp boundary between the 
low-density gas in the prompt axial outflow and the denser material expelled from the 
initial torus is very obvious. In all simulations presented here (and also others which 
probed the parameter range of the constant Bq), the MRI-driven accretion stream 
always overwhelms this early outflow. This seems inevitable, if only because the prompt 
axial outflow, by its very nature, is a transient effect. This transient is not to be 
confused with the funnel- wall and axial jets which are long-lived, stable features that 
arise once the turbulent accretion disk is established. 




2 4 6 8 10 12 



Figure 3: Interaction between nascent accretion flow from torus and axial outflow from the wind- 
up of the initial large-scale field. This image is produced from M u simulation data at t/M = 240. 
Gas density is plotted using a logarithmic colour scale. The arrows show the local direction of the 
poloidal component of the magnetic field. The figure spans a physical extent of 12 M; the black 
hole is centered on the left edge of the panel. 
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Once the accretion flow is established, the overall structure of the disk resembles 
previous results: the axial funnel, the funnel-wall jet, the extensive corona, and the 
main disk body are generic features of the simulations, with or without an initial large 
scale field. However, differences do arise in the corona and main disk due to the presence 
of the large-scale field and its interaction with the differentially rotating gas in the 
main disk body. Since the MRI operates within the initial torus, which is completely 
magnetized in these simulations (except for the Z simulation), the MRI has a much 
more extensive region in which to operate. The most striking manifestation of this is 
seen in plots of density, where extensive channels are observed; see Figure 01 These 
channels span a much greater radial extent in the M and S simulations than in the W 
simulation. 
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Figure 4: Comparison of extensive channels in the Weak-field (W; left panel) and Moderate-field 
(M u ) simulations at t/M = 1400. The presence of extensive channels at large radii in the main 
disk is a signature of the simulations with moderate to strong initial large-scale fields. Gas density 
is plotted using a logarithmic colour scale. The figures span a physical extent of 70 M; the black 
hole is centered on the left edge of the panel. 



3.3 The Strong-field Case 

In the S simulation, in which the large-scale field effectively washes out the initial 
poloidal loops, the evolution of the initial torus still proceeds along the lines of the 
other simulations, with similar late-time structures emerging. However, this simula- 
tion features much more violent, large-scale motions of gas in the disk and corona. 
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Animations show extended density channels feeding accreting gas onto the black hole. 
Those channels lying near the funnel wall can spontaneously break up into large vor- 
tices which move away from the black hole. Some of these vortices may then fall back 
and be accreted. These vortices appear to evolve, in part, due to the strong velocity 
shear along the funnel wall; they also show a vortical structure in the poloidal field, 
as shown in Figure [5J Though smaller versions of these vortices are observed in all 
simulations, the S simulation is by far the most spectacular. 
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Figure 5: Snapshots of density vortices moving away from the black hole following the break-up 
of a dense channel in the S simulation at times (left to right): t/M = 1800, 1808, 1816, and 1824. 
The arrows show the local direction of the poloidal component of the magnetic field in and around 
the vortices; the region surrounding the vortices has not only strong velocity shear (the outflow in 
the funnel is highly relativistic, whereas gas flow in the corona and funnel wall is much slower), 
but the poloidal component of the magnetic field also features a reversal on either side of the dense 
stream of gas containing the vortices. Each panel spans a physical extent of 15 M by 20 M; the 
black hole is located near the upper left corner in each image. 

Snapshots of the late-time state of the simulations suggest that the evolution of 
the S simulation is somehow retarded compared to the other simulations; by t/M = 
2000, it still exhibits extensive channels in the inner region whereas the W and M 
simulations exhibit a greater degree of turbulence. To investigate this, the S simulation 
was extended to t/M = 4000. Figure compares the density at t/M = 2000 and 
t/M = 4000. It is clear that the main disk eventually becomes turbulent, and that a 
bi-conical funnel and an extensive corona also emerge. This reinforces the notion that 
the outcome of the simulations is generic, but also suggests that the presence of the 
large-scale initial field alters the time scale of the evolution somewhat. By the time 
the S simulation is fully turbulent, the other simulations, especially the W simulation, 
are showing decreased MRI activity. 
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Figure 6: Comparison of density at t/M = 2000 (left panel) to t/M = 4000 (right panel) in the 
S simulation. Gas density is plotted using a logarithmic colour scale. The figures span a physical 
extent of 70 M; the black hole is centered on the left edge of the panel. 

4 Comparing the Simulations 

In order to compare and contrast the four simulations, records of shell-averaged scalars 
(e.g. density and pressure), as well as shell-averaged fluxes (e.g. mass, energy, magnetic 
stress) are available to provide a view of the simulations at high time-resolution (these 
quantities were computed at every 1 M increment in time during the evolution of each 
disk). In addition, time-averaging of code variables from the full data dumps provides 
a complementary view, by highlighting the more persistent features of the simulations. 



4.1 Mass and Energy 

The mass and energy fluxes ((p U r )(r, t) and (T r t )(r,t)) quantify the effect of the large- 
scale initial field on the accretion rate (at r m j n ), and the ejection rate of material at 
large radii. All integrated quantities are computed while distinguishing between bound 
{—hUt < 1) and unbound material; this simple distinction helps separate disk material 
(bound) from coronal/jet material (unbound). 

Table 13 shows the mass and energy passing through the inner boundary and ejected 
from the simulation volume through r/M = 100. For reference, the "Z" (zero large- 
scale field) simulation is also shown to provide a benchmark. 
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Table 2: Mass and Energy Accretion/Ejection. 



Simulation 


A Mhni,nH/Mn 


/\ 1V1h«1-.«iivi^ 

i,J -unDouna 


/M n 


/ \ r n i / xljn 


^unbound 


/En 


A ppyDTinn • 
.rVlst^l cllUil. 














Z (See Note) 


1.73 x 10 


2.84 x 10 


-4 


i r o ^ , 1 a— 2 

1.53 x 10 


o ot - , 1 r\- 

3.27 x 10 


-4 


T T T 


1.54 X 10 


3.34 x 10 


-4 


1 oc x^ i n — 2 

1.35 X 10 


3.93 x 10 


-4 




1 /i/i x ✓ in — 2 
1.44 X 10 


2.0/ x 10 


-4 


i op x^ i n — 2 

1.26 x 10 


a x^ 1 n- 

4.58 X 10 


-4 


M d 


o 11 w 1 n — 2 
2.11 X 11) 


o x^ 1 n- 

2.25 x 10 


-3 


1 to x^ 1 n — 2 

1.78 x 10 


o n*7 x^ i n- 

2.9/ x 10 


-3 


S 


2.66 x 1(T 2 


5 73 x 10 


-3 


2.25 x 10~ 2 


7.16 x 10 


-3 


Ejection: 














Z 


7.57 x 1(T 4 


1.68 x 10 


-3 


7.73 x 10" 4 


2.14 x 10 


-3 


w 


6.28 x 1(T 4 


1.37 x 10 


-3 


6.39 x 10" 4 


1.80 x 10 


-3 


M u 


5.02 x lO" 8 


3.31 x 10 


-4 


2.96 x HT 6 


1.81 x 10 


-3 


M d 


5.19 x 1(T 8 


4.83 x 10- 


-4 


3.10 x 10~ 6 


2.03 x 10 


-3 


S 


3.94 x 10~ 8 


8.13 x 10" 


-4 


3.75 x KT 6 


3.58 x 10- 


-3 



Note: This so-called zero-field simulation is a baseline simulation where the large-scale initial field is switched 
off, leaving only the initial poloidal loops with = 100. 



The mass/energy accretion and ejection numbers are obtained by integrating the 
corresponding shell-averaged fluxes over the full simulation (details in Appendix . 
and normalizing each to the initial mass/energy of the torus. The accreted mass and 
energy is seen to increase with the strength of the initial large-scale field, though there 
is some variability in the models. The rates for models Z, W, and M u tend to lie fairly 
close to one another. There is a noticeable jump in the accreted unbound mass and 
energy for models M d and S, with the total accreted mass (bound plus unbound) in S 
increased by a factor of ~ 2 over the Z, W, and M u simulations. 

The total amount of energy (bound plus unbound) ejected through r/M = 100 
is slightly enhanced by the initial large-scale field, while ejected mass is lower in the 
M and S simulations than it is in the Z and W cases. The nature of the ejected 
mass/energy, i.e. the portion that is bound vs unbound, shows a sharp shift in the M 
and S cases, indicating that more of the ejected material is truly ejected (has escape 
velocity) and would not return into the simulation volume should the outer radial 
boundary be pushed outward. 

Mass accretion as a function of time, (p U r )(r ra i n , t), is shown in Figure This 
figure departs from the convention adopted in this paper and plots simulation time 
in terms of orbital period at the initial pressure maximum. This highlights several 
important features arising from the initial large-scale field: 

• In all simulations, the accretion rate does not pick up until the first orbit is 
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complete: accretion is driven by the MRI operating within the disk. 

• The Z and W models show the characteristic spike at t = 1 orbit, marking the 
arrival of the dense accretion flow at the inner boundary. 

• The M and S models lack this sharp spike and show a more gradual increase, due 
in part to the "duel" between the early, tenuous outflow powered by the wind- 
ing up of the initial large-scale field outside the initial torus and the equatorial 
accretion stream from the torus. 

• The accretion rate of the Z and W models begins to die down after three orbits, 
indicating decreased MRI activity. 

• The accretion rate in the M u model shows some diminution after the third orbit, 
but remains stronger than the W and Z models. 

• The accretion rate in the and S models shows no sign of decreasing after 
the third orbit, and in fact seems to be growing in strength; the large spikes 
correspond to the arrival of dense channels and/or inbound vortices, similar to 
those shown in Figure 03 




No field 




Orbits 



Figure 7: Mass accretion rate at the inner radial boundary; each curve is plotted on the same 
scale, but shifted vertically for clarity (the dashed horizontal lines mark the zero for each curve). 
The solid curves mark the total accretion rates (bound plus unbound); while the dashed curves 
show unbound material only. Most of the accreted material is bound. 
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4.2 Time-averaged Data 



The overall structure of the main disk, corona, and funnel is best viewed in time- 
averaged plots, which smooth out the often intense variability of the simulations. The 
three main differences between the simulations are seen in the density distribution in 
the corona and axial funnel, in the overall pressure distribution, and in the relative 
strength of the toroidal and poloidal components of the magnetic field. 

Figure |H1 shows the density distribution in each simulation, averaged from t/M = 
1200 to 1600. This figure shows the more expansive region occupied by the main disk 
in the M and S simulations, due to the action of the MRI on the more extensive initial 
magnetic field. The figure also shows that the M and S simulations tend to produce 
a narrower funnel, and hence a more collimated axial jet; the figure suggests that a 
tight, nearly cylindrical outflow is achieved. However, Figure HO suggests that this is a 
transient phenomenon: the magnetic structure at large radii appears to evolve over a 
longer time in the M and S simulations, and the tight collimation eventually vanishes. 



-10 -8 -6 -4 -2 




Figure 8: Comparison of time-averaged (from t/M =1200 to 1600) density for all simulations, from 
left to right: W, M u , M<j, S. Logarithmic colour scale is normalized to initial maximum density. 
Figures span a width 50 M and a height of 100 M; the black hole is centered on the left edge of the 
panel. 
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Figure IH1 shows the pressure distribution in each simulation, averaged from t/M = 
1200 to 1600. Both the total pressure (gas plus magnetic) as well as the magnetic 
pressure are shown. The W simulation has a total pressure distribution essentially 
identical to simulations with no initial large-scale magnetic field, namely a spherical 
"halo" structure in the funnel and a horizontally stratified structure in the main disk. 
The M and S simulations show a marked departure from this, with a significantly 
stronger magnetic pressure in the corona and at large radii in the funnel. The region 
of high magnetic pressure in the coronal region immediately above and below the inner 
torus in the and S simulations has a scale height greater than that of the main 
disk. Evidence that the M and S simulations are less evolved at large radii is also 
found in these plots: the pressure distribution at large radii still retains some of its 
initial structure. A region of high magnetic pressure along the spin axis is likely an 
artifact of the reflecting boundary conditions used in axisymmetry. 




Figure 9: Comparison of time-averaged pressure (from t/M =1200 to 1600) for all simulations, from 
left to right: W, M u , Mj, S. Top row shows magnetic pressure; bottom row shows total pressure 
(gas and magnetic). The logarithmic colour scale is normalized to maximum total pressure at 
t/M = 1600 for simulation S. The figures span a width 50 M and a height of 100 M; the black hole 
is centered on the left edge of the panel. 
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4.3 Evolution of the Magnetic Field 



Since the material in the corona and the axial jets is entirely unbound (—h Ut > 1), the 
diagnostic tracking the evolution of magnetic pressure in unbound gas provides a view 
of the time-evolution of the magnetic field in these regions. Figure ITU1 shows the time- 
evolution, for all four simulations, of unbound magnetic pressure on the shell passing 
through the initial pressure maximum, a region where Figure El indicates a pronounced 
difference between simulations. The M and S simulations show a prompt rise in 
unbound magnetic pressure, due to the stronger large-scale initial field, which is wound 
up by frame-dragging. The growth of unbound magnetic pressure in the W simulation 
lags considerably, reinforcing the idea that this simulation is driven by the MRI acting 
on the initial poloidal loops, that is that unbound material is produced during the 
saturation phase of the MRI, as it is in "conventional" disk simulations. Furthermore, 
after the initial sharp rise, the magnetic pressure for the M and S simulations continues 
to increase steadily with time, whereas the magnetic pressure in the W simulation 
reaches a peak at t/M ~ 850, and then decreases. 




Figure 10: Time-evolution of the shell-averaged, unbound magnetic pressure at r/M = 16.1. 
Unbound material is found in the corona and the axial jets, and the radial distance of the shell 
shown here passes through the initial pressure maximum. The solid line denotes the W simulation; 
the dotted line denotes M u ; the dashed line, Ma; the dash-dotted line, S. 
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In the M and S simulations, there is a sharp jump in magnetic pressure at t/M ~ 
300, with the greatest jump occurring for the S simulation. A snapshot of the magnetic 
field at this early time is given in Figure ITT1 This figure shows the toroidal component 
of the magnetic field plotted in grey scale, along with arrows showing the local direction 
of the poloidal field. The S simulation stands out in this figure: the region between the 
thin equatorial inflow and the funnel wall is occupied by a region of intense, smooth 
toroidal field. In addition, while the W and M models show regions of strong radial 
streams of toroidal field near the initial pressure maximum, indicating the growth phase 
of the MRI, the S simulation shows exactly the opposite, with a region of weaker radial 
streams at much greater radii. The sequence of images suggests that the stronger initial 
large-scale field generates a sheath of toroidal field around the torus which squeezes 
the region of MRI activity to larger radii. This reinforces the notion that the evolution 
of the S simulation is in many respects retarded compared to the other simulations. 
Another important feature in this figure is the presence of toroidal field in the axial 
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Figure 11: Snapshot of magnetic field at t/M = 330 for all simulations, from left to right: W, 
M u , Mj, S. The arrows show the local direction of the poloidal component of the magnetic field; 
the toroidal component is shown in grey-scale, in code units. Each panel spans a physical extent 
of 30 M; the black hole is centered on the left edge of each panel. 

funnel, with the strongest funnel field found in the S simulation. Animations show 
that this toroidal component of the funnel field is highly time- variable, appearing in 
bursts that suggest the launching of torsional waves through the funnel. These bursts 
originate near the inner radial boundary, indicating an origin in the ergosphere. 

Finally, to discuss the onset of turbulence in the main disk, take the b r component 
of the shell-averaged angular momentum flux (bound component), ((b r b^) \-hUt>l)(T, t). 
Figure IT2l shows spacetime diagrams for this quantity for each of the four simulations. 
The W model shows a steady growth in the first few 100 M of simulation time, followed 
by a quasi-steady state, then a decrease after 1500 M. The region of elevated magnetic 
stress is confined to the vicinity of the initial pressure maximum, where the initial 
magnetic field is concentrated. The decrease in stress at late times is consistent with 
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the decrease in accretion rate seen in Figure [7] suggesting a decline in MRI activity. 
The M and S simulations show two important differences from the W simulation: the 
onset of elevated stress is delayed, with the delay most pronounced in the S simulation; 
and the region of elevated stress expands with time to large radii, engulfing a greater 
portion of the disk. The region of elevated stress extends to the end of the simulation 
in the and S simulations, suggesting that the turbulent phase of the MRI is much 
longer-lived when a strong initial large-scale field is present. 

-0.0030 -0.0020 -0.0010 0.0000 0.0010 0.0020 0.0030 




Figure 12: Spacetime diagrams of magnetic stress (bound material), ((b r b^) \-hUt>i)( r i *)> from 
left to right: simulations W, M u , M^, S. 

The duration of the growth phase for stress in the S simulation can investigated by 
considering the records from the extended S simulation. Figure fTTA shows the extended 
diagnostic for this simulation: the period of growing/elevated stress, though consid- 
erably longer in the S simulation, is nonetheless finite in duration. This simulation 
shows a decrease in intensity after t/M ~ 2300. 
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Figure 13: Spacetime diagram of magnetic stress (bound material), ((b r b^) |-/i£/ t >i )(r, t), for 
extended S simulation. 

5 Discussion 

This paper has presented results of a group of axisymmetric GRMHD simulations of 
accretion disks embedded in a large-scale magnetic field. To investigate the role of this 
large-scale field, all simulations used a common initial state: Kerr black hole with spin 
a/M = 0.9, equilibrium torus perturbed by small, random density fluctuations, and an 
MRI seed field consisting of weak poloidal loops confined to the core of the torus. The 
simulations differed only in the scaling constant for the initial large-scale field. This 
initial field was termed weak, moderate, or strong, in relation to its effect on the loops 
of poloidal field: weak if it left the loops intact; moderate if it cancelled the portion of 
the loops anti-parallel to it, on either side of the initial pressure maximum; and strong, 
if it overwhelmed the loops. 

The large-scale initial magnetic field does not alter the general features of the late- 
time structure of the accretion system: the features identified in DHK (main disk, 
inner torus, plunging region, coronal and axial funnel) remain qualitatively the same. 
However, the relative strength of the initial large-scale field was found to have an effect 
on the details of the evolution of the initial torus. The key points are the following: 

• the stronger the field, the greater the accretion rate, though the effect is modest; 

• the stronger the field, the greater the ejection rate of energy, though the effect is, 
again, modest; 

• the stronger the field, the greater the fraction of ejected material that is unbound; 
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• the stronger the field, the greater the scale height of the region of elevated mag- 
netic pressure in the corona; 

• the stronger the field, the later the onset of the growth phase of the MRI; 

• the stronger the field, the more intense the density channels prior to the onset of 
turbulence. 

In addition, the large-scale field, due to the higher coronal magnetic pressure, tends to 
produce a more collimated jet, at least in the early stages, effectively by squeezing low- 
density funnel gas closer to the spin axis (Figure ITT|) . This suggests that a cylindrically- 
collimated jet is produced very close to the black hole, though simulations with a larger 
radial grid should be carried out to verify whether this collimation is found at large 
distances. Indications here are that the collimation is confined to the region near the 
black hole, and that it is short-lived; a more conical shape is achieved at late times in 
the extended S simulation. 

The mass and energy accretion and ejection rates in Table [21 can be compared to 
those found in earlier simulations, notably those in DHK and D05a. It is possible to 
compare the Z simulation to its closest D05a counterpart, the 3D model KDPi (KDP 
with r = 4/3). The rates are found to be systematically lower in the simulations 
discussed here, by a factor of ~ 2. The most important difference which can account 
for the lower rates is the use of axisymmetry, which tends to reduce the accretion rate 
(De Villiers & Hawley, 2003; DH03b). The rates given in TableEJcan also be compared 
to the 3D simulation in D05b, model -R v f 3D ; which is very similar to the W simulation 
discussed here. The mass accreted onto the black hole for this earlier simulation is 
lower than that found here, by a factor of ~ 5, whereas the ejected quantities are 
systematically higher than those shown here, again by a factor of ~ 5. 

Hawley & Balbus (1992) discussed the evolution of a shearing sheet containing an 
initially weak vertical magnetic field, and documented the emergence of the channel 
solution, radially elongated regions of elevated density sandwiched between regions of 
low density and strong magnetic field, as a direct consequence of the MRI interacting 
with the vertical field. In subsequent studies of axisymmetric disks, these channels have 
emerged as a characteristic feature (channels also arise in 3D, but quickly dissipate; 
DH03b). Channels are found in the set of simulations discussed here, and the channels 
are enhanced in strength and radial extent by the presence of the pervasive large-scale 
field. This is not surprising: the large-scale initial field permeates the entire initial 
torus, and the MRI is thus provided with a much larger initial region in which to 
amplify magnetic field. The simulations with a stronger initial large-scale field give 
rise to a more vigorous evolution, a feature better appreciated by viewing animations 
of gas density, though the figures presented here capture some of this detail. 

What may seem paradoxical is that the strongest initial field simulated here shows 
the latest onset of turbulence: though the S simulation is by far the most intense, the 
growth and subsequent saturation of the MRI are the slowest of all the simulations. 
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However, this is not a new finding: HB92 (c.f. Fig. 5) observed the same effect in 
the shearing sheet simulations, namely that stronger initial fields show a delayed onset 
of the growth phase of the MRI when all other simulation parameters are kept con- 
stant. HB92 also found that stronger initial magnetic fields generate a much stronger 
maximum field, provided that the initial field was not so strong as to be stable against 
the MRI. The current simulations are consistent with these earlier findings, for the 
same reasons: though dealing with a more complicated initial state, the Z, W, M, and 
S simulations differ only in the strength of the initial large-scale field, and all other 
parameters (grid size, time step size, disk configuration) are identical. In addition, 
the strength of the initial large-scale field was kept to a reasonable level, that is with 
f3 3> 1 in the initial torus, in recognition of the fact that the importance of the MRI 
as the mechanism of angular momentum transport in astrophysical disks rests on the 
fact that it operates on weak magnetic fields. Simulations that begin with inordinately 
strong magnetic fields may well preclude the development of the MRI and may, as was 
observed in some early tests for the present study, severely disrupt the initial torus to 
the point where the physical plausibility of the simulation is in doubt. 

Another interesting finding is the noticeable difference in outcomes between the 
and M u simulations which differ only in the relative orientation of the large-scale initial 
field. The outcome of the M u simulation, where the field is in the usual "up" sense, is 
very similar to the weak-field W simulation, whereas the simulation, where the field 
is oriented "downward" , bears a closer similarity to the strong-field S simulation. The 
different outcomes seem to be rooted in the interaction between the large-scale initial 
field and the poloidal loops: the M u simulation disrupts the poloidal loops in such a 
way as to strengthen the field in the narrow region near the inner edge of the torus, 
but it also significantly weakens the field at the center of the torus (just outside the 
initial pressure maximum). In this sense, the M u simulation begins from a weaker net 
initial field. The simulation does exactly the opposite, by reinforcing the poloidal 
loops at the center of the torus it produces a relatively stronger initial field. 

There are many aspects to these simulations that were only briefly touched upon in 
this survey. More detailed analysis of the simulation data should be undertaken, as well 
as complementary simulations to probe the effect of other simulation parameters, most 
importantly black hole spin. Further investigation is warranted in the properties of the 
funnel: the transient collimation near the base of the funnel due to enhanced coronal 
pressure; the appearance of bursts of toroidal field in the funnel and their presumed 
association to torsional Alfven waves; and the behaviour of the axial jets at large radii. 
More importantly, the lifting of the axisymmetry constraint by proceeding to full 3D 
simulations is essential to answer questions regarding the long-term properties of the 
accretion disks when immersed in large-scale fields. 
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A Equations of GRMHD 



The GRMHD code, described in its original form in DH03a, uses a finite difference 
approach to solve the equations of ideal general relativistic magneto-hydrodynamics in 
the spacetime of a Kerr (rotating) black hole in Boyer-Lindquist coordinates, (t, r, 6, (f>). 
The GRMHD code solves the equation of continuity, \/^(pU^) = 0, the energy- 
momentum conservation law, V M T^ v = 0, and the Maxwell's equations, V fl i ? ' 11 ' = 
4-7T J v and V^*-F^ = 0, which, in the MHD approximation, reduce to the induc- 
tion equation d$ F a p + d a Fp$ + dp Fs a = 0. In the above expressions, p is the 
density, the 4- velocity, T^ u the energy-momentum tensor, and F a p the electro- 
magnetic field strength tensor. The energy momentum tensor is given by T^ u = 
(ph + ll&f) U» U u + {p + ) - W b v where h = 1 + e + P/p is the specific 
enthalpy, with e the specific internal energy and P = pe(T — 1) the ideal gas pres- 
sure (r is the adiabatic exponent); ||6|| 2 = b^b^ is the magnetic field intensity; and 
b^ = *F iiV U v /(A-k) is the magnetic field 4- vector. The induction equation is rewritten 
in terms of the Constrained Transport (CT; Evans & Hawley, 1988) magnetic field vari- 
ables & = e ijk F jk , as dtB l -dj [V 1 B j - V j B 1 ) = 0, where V 1 = W /U* is the transport 
velocity, and U* = W/a, with W the Lorentz factor. The CT algorithm ensures that the 
constraint diB 1 = is satisfied to rounding error. The Kerr metric is expressed in Boyer- 
Lindquist coordinates, for which ds 2 = gu dt 2 + 2 g t< ^ dt deft + g rr dr 2 + ggg d9 2 + g^ d(f) 2 ; 
a = (—g u ) is the lapse function. Geometrodynamic units are used, where G = c = 
1, and time and distance are given in units of the black hole mass, M. 

A recent effort at development and testing of the GRMHD code has yielded im- 
proved stability in highly magnetized flows by adding artificial diffusion terms to the 
equations of continuity, energy conservation, and momentum conservation. Though 
a detailed discussion of these tests is beyond the scope of this paper, the reason for 
adding such diffusion terms is quite simple: truncation error due to finite differencing 
of time derivatives in expressions of the form 

dtX + ^pdji^V'X) (1) 
v T 
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gives rise to grid noise that may be amplified to unacceptable levels over time. The 
remedy is to add compensatory flux terms which cancel the leading term in the remain- 
der of a formal Taylor expansion of the time derivative operator. The time derivative 
term, expanded to second order, is given here in schematic form: 

yn+l -yn 

' At 1 ~f (2) 

where the superscripts denote time levels and the subscript denotes a spatial grid index. 
Using the second-order term can be rewritten as 

d l X ^l_ dj [ V i{d i ^jV i X)} (3) 
and fed back in to the discretized form as an artificial diffusion term: 



At ~ j 



-3 



(4) 



Numerical tests have shown that the diffusion constant a = 6 x 10 yields good results. 

These modifications were important for the simulations discussed here since the 
addition of the initial large-scale magnetic field produces conditions near the black 
hole where shocks in a highly magnetized flow can give rise to numerical artefacts. 
However, it was also found that these modifications to the GRMHD code had only a 
minor impact on "conventional" accretion disk simulations (i.e. those without an initial 
large-scale magnetic field) such as the KD models discussed in DHK and companion 
papers. This is corroborated by the agreement mentioned in ^between the present 
simulations and those in D05a. The version of the GRMHD code used to produce 
DH05b, an intermediate form of the code, tended to overestimate funnel energetics. 



B Details on the Initial State 

The initial torus, discussed in detail in DHK, is an equilibrium disk solution to the 
equations of GR hydrodynamics which is then pertubed with a small, random density 
fluctuation and overlaid with a magnetic field that serves as a seed field for the MRI. 

The equilibrium solution is found by considering a disk with a power-law rotation, 

n = 7]\- q (5) 

where 

i] is a constant, and q a positive number. The hydrodynamic momentum equation is 
(Hawley, Smarr, & Wilson, 1984) 



dt(S j ) + ±=:d i ^{S j V i ) + l - d j gf" + ad j {P)=0, 



(7) 
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where Sj = phW Uj is the momentum. 



The momentum equation is simplified by imposing time- independence, axisymme- 
try, and requiring that there be no poloidal motion. This simplified equation reads, 

<*d 3 (P) + \ d j9 ^ = 0. (8) 

Use the definition of the momentum 4-vector and the definition of specific angular 
momentum, I = —U^/Ut, to obtain 



d< (P) U t 2 



d 3 {U t - 2 ) + U t 2 (-d 3 - g** + l d, g^) d 3 I, (9) 



ph 2 
where U t ~ 2 = g u - 2 1 g t( t> + I 2 g**. 

Assume constant entropy, Tds = 0, so that dh = dp/ p. Use © and (jHJ) to write 
I = n\ 2 = c\ 2 - q and Q = tr 2 ^ -2 ) f«/(«~ 2 ) = kl a . So © can be written in integral 
form 

'■ h dh i r Ut d(ur 2 ) r l ki a 



i > i , , ^dl (10) 

h h 2Ju. Ur 2 J i l-kl^ 1 v ' 

where a = q/{q — 2) and Xi n refers to the quantity in question evaluated on the surface 
of the disk. Clearly, hi n = 0, and a general solution admitting a choice of surface 
binding energy Ui n is given by 

h( r 9) = Uin f( lin ) (11) 

u t (r,e)f(i(r,e)y (iL) 

where f(l) = \\l - kl a+1 \\ 1/(a+1} or /(fi) = ||1 - k~ l l a ft(«+i)A*|| 1/(a+1) . Using the 
equation of state and the definition of enthalpy, the internal energy of the disk is 

r n\ If Ui n f(li n ) \ /io\ 

^'f UwiM) -')' (12) 

For a constant entropy adiabatic gas the pressure is given by P = pe (T — 1) = K p r , 
and density is given by p = [e(T — I) / K] l ^ T ~ l \ 

For the Schwarzschild metric, these analytic relations completely specify the initial 
equilibrium torus. In the Kerr metric, A 2 is an implicit function of I, requiring an iter- 
ative approach to solve for the equilibrium state. However, by using the Schwarzschild 
expression for A 2 in the Kerr metric, the above results can be used to initialize a disk 
that maintains a sufficiently good equilibrium to use an initial state for numerical 
studies of the MRI. 

Poloidal loops of magnetic field are laid down along isodensity surfaces within the 
torus by defining = (^,0,0,^), where 

a - J k (P ~ Pcut) f or p > pcut n on 

^ (1 ° OPS) -\0 for p<p cut > (13) 



25 



where p cut is a cutoff density corresponding to a particular isodensity surface within 
the torus. The constant k is set by the input parameter /3, the ratio of the gas pressure 
to the magnetic pressure. The constant p cu t = 0.25p max is chosen to keep the initial 
magnetic field away from the outer edge of the disk. 

The large-scale magnetic field is introduced in a very straightforward manner by 
adding Wald's solution (W74) to the poloidal loops: 

A 4> = A 4> (loops) + -y (9<M> + 2ag t(j) ) (14) 
The magnetic field is initialized using B r = —dgA^ and B 6 = d r A^. 



C Evolution Diagnostics 



The code diagnostics consist of shell- averaged and volume-integrated quantities com- 
puted at regular time intervals. The average of a quantity X on a shell, radius r, 

(X)(r,t) = [ [x^=gd0d<f> (15) 



where the area of a shell is A(r) and the bounds of integration range over the 9 and 
(f> grids. For these simulations, shell-averaged values of density, (p), angular momen- 
tum, (pi), gas pressure, (P), magnetic pressure (||6|| 2 )/2, angular momentum (phU^) 
and binding energy (phUt) are computed. Two sets of shell averages are available, 
based on whether the fluid at a given point on the shell is bound (—hUt < 1) or un- 
bound; for instance, the history of the unbound component of quantity X is written 
(X\_ hUt>1 )(r,t). 

Fluxes through the shell are computed in the same manner, but are not normalized 
with the area. The rest mass flux (pU r ), energy flux 

(T r t )(r,t) = (phU r U t )(r,t) + (\\b\\ 2 U r U t )(r,t) - (b r b t )(r,t), (16) 

and the angular momentum flux 

CZ%)(r,t) = (phirUJfrt) + <||6|| 2 Cr^)(r,t) - <& r ^)(r,t) (17) 

are computed. To ensure flexibility, each of the three components in the above sums 
is stored in a separate file. In addition, bound and unbound components of fluxes are 
computed and stored separately. 

Volume-integrated quantities are computed using 

[Q](t)=j j j Q^g-drd6d4>. (18) 

The volume-integrated quantities are the total rest mass, [pU l ~\ , angular momentum 
[r*^] , and total energy [T* t ] , distinguishing between bound and unbound material. 
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